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Abstract 



A new hybrid scheme for numerical relativity will be presented. The scheme will employ a 3-dimensional 
spacelike lattice to record the 3- metric while using the standard 3+1 ADM equations to evolve the lattice. 
Each time step will involve three basic steps. First, the coordinate quantities such as the Riemann and 
extrinsic curvatures are extracted from the lattice. Second, the 3+1 ADM equations are used to evolve the 
coordinate data, and finally, the coordinate data is used to update the scalar data on the lattice (such as 
the leg lengths). The scheme will be presented only for the case of vacuum spacetime though there is no 
reason why it could not be extended to non-vacuum spacetimes. The scheme allows any choice for the lapse 
function and shift vectors. An example for the Kasner T 3 cosmology will be presented and it will be shown 
that the method has, for this simple example, zero discretisation error. 

1. Introduction 

In a recent paper [1] a new numerical scheme was developed and successfully applied in the 
construction of time symmetric initial data for a Schwarzschild black hole. The method was 
based on two key ideas, one to use short geodesic segments as the legs of the lattice and two, 
to employ Riemann normal coordinates local to each vertex. 

Our primary aim in this paper is to present a hybrid scheme in which a 3-dimensional 
spacelike lattice can be evolved using the standard ADM vacuum 3+1 equations. 

The general scheme envisaged in [1] can be summarised as follows. The lattice is a network of 
vertices and legs and is fully specified by the connectivity matrix and the leg lengths. To each 
vertex there is an associated computational cell defined by some local sub-set of the lattice 
(eg. the legs and vertices of the tetrahedra attached to this vertex, though larger subsets 
could be considered). Riemann normal coordinates are employed in each computational cell 
leading to a metric of the form 



where is chosen to be diag(l, 1, 1) and e is a typical length scale for the computational 
cell. The crucial assertion that the legs are geodesic segments of this metric leads directly 





to the following equation for the leg length Ly between vertices (i) and (j) 

L% = g^Ax^Ax^ - x»x\x a jx] + 0(e 5 ) (1.2) 

where Ax^ = Xj — xf. For each computational cell there would be a small set of such 
equations, which in [1] where referred to as the smooth lattice equations. These could, in 
principle, be solved as a coupled set of equations for (estimates of) the Riemann tensor and 
any coordinates not fixed by gauge freedoms. The computations for any one cell are fully 
decoupled from those of any other cell. Thus this is a local problem and should be amenable 
to Newton-Raphson methods. It is at this point that the design of the lattice becomes 
crucial. It must be such that the equations in each cell yield unambiguous estimates for 
the Riemann tensor. If the system is singular then there will exist continuous families of 
solutions for the curvatures. This indicates that the lattice lacks sufficient information to 
tie down a unique estimate (or at worst a finite discrete set of estimates) for the curvatures. 
The solution then would be to add extra information to the lattice such as extra leg lengths. 
The lattice chosen in [1] did not suffer from this problem. 

However, subsequent experiments [2] have shown that many seemingly suitable lattices 
(such as in Fig 4) are singular. Following these experiments a slight variation to the above 
scheme has proved to be very successful. It entails the addition of extra information in the 
form of the polar directions of the legs attached to the central vertex of the cell. Since the 
connection vanishes at the origin (by definition of the RNC) this is equivalent to specifying 
the spherical polar angles of the vertices at the ends of these legs. This not only supplies 
the missing information but it also fully decouples the set of smooth lattice equations. 

The modified scheme works as follows. Each computational cell will, for the remainder of 
this paper, be assumed to consist of the central vertex and the tetrahedra attached to this 
vertex. For each cell there will be one central vertex, for which we will always use the label 
0, and a set of vertices on the surface of the cell. To each surface vertex (i) we assign the 
polar coordinates (#j,0j). With this, the equations (1.2) for the radial legs are satisfied by 

Xi = L i cos fa sin Q{ 

yi = Loi sin fa sin 9{ (1.3) 

%i — cos 6i 

where L i is the length of the radial leg from the origin of the cell to the surface vertex (i). 
The ( Xi,yi,Zi) are the Riemann normal coordinates for vertex (i). Next, the curvatures can 
be computed by solving, by a least squares method, the leg length equations (1.2) applied 
to the legs on the surface of the computational cell. Least squares, or some other method, 
will be necessary since these equations are almost certainly overdetermined. 
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In this scheme the basic data given on the lattice are the L{j,9i and fa from which we 
compute, by the above equations, the x^ and the R^ a f3- 

Note that for each leg on the surface of the cell there is an associated angle subtended at 
the origin of the cell. In terms of the spherical polar coordinates for the vertices (i) and (J) 
of the leg (ij), this angle, cty, can be computed from 

cos aij = cos 9{ cos 0j + cos(fa — <pj) sin 9{ sin 9j (1.4) 

and in terms of the Riemann normal coordinates from 

^Oi^oj r-\ 

cosaij = gyw— — (1.5) 

Both of these equations will be used later in this paper. 

Note that in the construction of our Riemann normal coordinates we have taken advantage 
of some coordinate freedoms. For example we have used translational freedoms to tie the 
origin of the RNC to the central vertex. We have also chosen our coordinate axes to be 
orthogonal at the origin (hence g^ v = diag(l, 1, 1)). However, there are some residual gauge 
freedoms, namely the three degrees of freedom to rotate the coordinate axes about the origin. 
This freedom can be accounted for by specifying three of the polar angles 9{ and fa, such as 
6>i, fa and #2- 

2. Smooth lattice kinematics 

In the familiar 3+1 formulation of general relativity a typical spacetime can be represented 
in two equivalent ways, one as an evolving 3-metric on 3-dimensional manifold and two, as 
sequence of distinct 3-dimensional hypersurfaces of the spacetime. In the later picture, each 
hypersurface records one instant in the evolution of the 3-metric. The information that one 
must specify to fully define the 4-metric on the spacetime are the 3-metric, the lapse function 
and the shift vectors on each hypersurface. In this section we shall adapt this picture to one 
in which each hypersurface is a smooth 3-dimensional lattice. 

The general idea will be to allow the basic data for each lattice, namely the L{j,9\ and fa, 
to be functions of time. On each hypersurface we will continue to employ Riemann normal 
coordinates and the above equations for computing the curvature components. Since the 
lattice data are now presumed to be functions of time so too must the coordinates xf and 
the curvature components R^ va i3- We will also impose our previously stated gauge conditions 
on each hypersurface, in particular that g^ v = diag(l, 1, 1) and that the origin of the RNC 
remains forever tied to the central vertex of the cell. One could choose a dynamical set of 
gauge conditions but to do so at this stage would be an un-necessary complication. 
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It is customary to use Latin indices to denote spatial components (eg. gij for the 3-metric) 
and Greek indices for spacetime quantities (eg. g^ v for the 4- metric). As we would like to 
use Latin indices to denote simplices of the lattice we shall choose to use Greek indices for 
spatial components. Thus in the following g^ u , R^ va i3 etc. will denote 3-tensors that live in 
each 3-dimensional hypersurface. 

We will be making frequent reference to the current and future hypersurfaces and the various 
structures that define and link these hypersurfaces. We will use So to represent the current 
hypersurface and T,g t to represent the hypersurface just slightly to the future of So- The 
time coordinate t is constant on each hypersurface with t = to on and t = to + St on E^. 
The lapse function will be denoted by iV and the shift vector by N^. Both the lapse function 
and shift vectors will be defined on the vertices of the lattice. Note that even though we 
are using Riemann normal coordinates on each hypersurface the 4-dimensional coordinates 
(t, x^) will not, in general, be in Riemann normal form. 

Consider any 3-dimensional smooth lattice and focus attention on one vertex. Through this 
vertex one can construct three well defined curves - the timelike worldline of the vertex, the 
timelike geodesic tangent to the normal vector at the vertex and the timelike worldline of 
the observer with constant spatial coordinates. It is important to realise that these may be 
three distinct curves. Thus upon extending these curves from their common origin in So to 
T,st we can expect to obtain three new points in E^. This situation is depicted in Fig (1). 
In the usual 3+1 picture only two of these curves, the normal geodesic and the observers 
worldline, enter into the discussion. In the smooth lattice we must also consider the freedom 
for vertices to drift relative to these two curves. For this reason we will introduce a new 
vector, the drift vector 7 M . 

From Fig (1) we can clearly see that the shift and drift vectors at vertex (i) are related by 

dxf 



dt 



*- = - 7 f + iVf (2.1) 



where xf(t) are the time dependent coordinates of vertex (i). The dxf/dt can be easily 
computed from (1.3). Note that in our chosen gauge xf = dx^ /dt = at the origin of the 
RNC. 

One now has a choice as to which of the drift or shift vectors should be freely specified on the 
lattice. Since the drift vector has a strong geometric appeal we shall choose to freely specify 
it and to use the above equation (2.1) to compute the shift vector. A convenient choice is to 
set the drift vector to be zero at the origin (ie. the vertex at the origin of the RNC in Eo is 
evolved along the timelike normal). Though, of course, there may be occasions where such 
a choice is not appropriate. 
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We will now introduce a second lattice which will help us later on in making the transition 
to the standard ADM 3+1 equations. Consider one hypersurface and its associated lattice. 
Consider now a second lattice coincident with the original lattice. We will refer to these as the 
primary and shadow lattices respectively. The shadow lattice will be chosen so such that its 
vertices are evolved along the normals to the hypersurface (ie. as if the drift vectors were set to 
zero everywhere). We do this because the ADM 3+1 equations take on a particularly simple 
form when expressed as time derivatives along the normals. The shadow lattice will only be 
employed for one time step. Upon completion of the time step the current shadow lattice 
will be discarded and a new shadow lattice created coinciding with the updated primary 
lattice. Each of these shadow lattices are introduced only as an aid in the exposition of our 
algorithm - they need never be created in any computer program. 

Our current task is to see how we should specify the data on the shadow lattice so that it 
evolves along the normals. 

On the primary lattice we have various quantities such as Lij, 9i, (pi and xf all of which we 
assume to be functions of time t. Their counterparts on the shadow lattice will be denoted 
by the addition of a dash (ie. L^-^xj 1 etc.). The dash on the x'f 1 should not be confused 
with coordinate transformation - the x'f 1 are just the x^ coordinates of vertex [i'). On the 
initial hypersurface So all of the corresponding dashed and un-dashed quantities are equal 
(ie. L\j = Lijjx'f 1 = xf) though they may drift apart between successive hyp ersurf aces. 

By inspection of Fig (1) one can easily verify that 

dt dt h y } 



Consider now equation (1.2) applied to both lattices, 



Lij(t) — g^Ax^-Ax^ — -R^avp xfx'^x'jXj 



3 



All of the terms on the right hand side should be considered to be functions of time (except 
the g^ v which we choose to be diag(l, 1, 1)). Differentiating with respect to t and using the 
above equation for dx'f'/dt we obtain on So 

^ = - z^A 7 gA4 + W<* a A + (2-3) 

where A 7 g = - 7 f . 
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A similar argument can be applied to the angles a y -. Starting with 

Ax oi Ax oi 

cos ay (f) = g^—— — J - (2.4) 

L J oi*-'oj 

cosa^(t) =^ (2.5) 

L oi L oj 

we obtain, by direct differentiation, that on So 

da'jj _ dajj , cosajj / 1 d , 1 d , 

dt dt sin ay \L idt L j dt J 



LoiL j sin oiij 



{^oA<3 + ^o^oi) 



This can be simplified by using the above equation for dL' /dt and noting that in our gauge 
Xo = on Eq. The result is that on So 



da ij _ dajj 
~df ~ ~df 



(2.6) 



The dL'ij/dt and the da'^/dt measure rates of change along the normals to So. They should 
therefore be expressible in terms of the extrinsic curvatures K'^. 

Consider the shadow lattice and the coordinates x^. Let us introduce a second coordinate 
system x"^ such that the x"^ are constant along each worldline of the vertices of the shadow 
lattice. This clearly establishes a time dependent coordinate transformation between x^ and 
x"^. Note that on So both coordinate systems are identical. Note also that on subsequent 
hypersurfaces the x"^ need not be in Riemann normal coordinate form. 

Since the worldlines of the shadow lattice are tangent to the normals to So we must have 

^jlV _ O A T TSH 

~df ~ ~ 2NK ^ 

In the x"^ coordinates the L'^{t) and a^-(i) may be estimated from 

L%{t)=g'/ v {t)Ax^Ax'(J + 0{^) (2.7) 
L'oiL'oj cosak(f) = g'^Ax'^Ax 1 /," + 0(e 4 ) (2.8) 

where g'/ v is evaluated on the geodesic tangent to the normal to So at the origin of the RNC 
Since the x"^ need not be in Riemann normal coordinate form we can expect these estimates 
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to be less accurate than their Riemann normal coordinate counterparts (eg. equations (1.2) 
and (1.5)). This is the price we pay for not using Riemann normal coordinates. Increased 
accuracy could be obtained for the Lf-{t) by evaluating the right hand side at the centre 
of the legs. This then forces us to explicitly account for variations of K^ v throughout the 
cell. It is not clear how one might make a similar 'centred' approximation for the angles 
a'ijit). Rather than trying to resolve these issues in this paper we shall settle on the simple 
estimates given above. 

Differentiating the previous two equations with respect to time and using the above equation 
for K'l^ we obtain, on So, 

dL' 

L'ij-jf = -(NK^oAxf^j (2.9) 

j t {L' 0l L' 0j cos ay = -2(NK^) Ax» oi Ax» oj (2.10) 

This system of equations for the six K^ v at the origin of the RNC will, for almost all 
lattices, be highly overdetermined (there being many more Ly and cty than the six K^). 
This problem can be overcome in one of two ways. We could reject all but six of the above 
equations. Or we could use all of the equations in a least squares estimation of the six K^ v . 
We will assume in the following that the least squares method has been used. In fact we 
will find that this issue of overdeterminism will crop up on a number of occasions in this 
paper. In all cases we will resort to a least squares solution as it is systematic and easy to 
implement in a computer program. 

At this point we are in a position to calculate both the extrinsic and intrinsic curvatures 
tensors, at the origin of the RNC, given just the basic lattice data, namely the Lij, 6i, 
their time derivatives and for any choice of chosen lapse function and shift and drift vectors. 
However we have not discussed how K^ v \ a and may be calculated (both of which will 
be needed in the next section). 

Consider first the calculation of N^ u . Since N is a scalar function we have = dN/dx^ 
which we can approximate at the centre of the leg (ij) by 



= (2.11) 



Ni — N4 

Such an expression can be formed at the centre of each leg of the computational cell. From 
this data we can form a linear approximation N^(x a ) to N^(x a ) 

N^{x a )=N^ + N^x» (2.12) 

where the constants and N^ v = N\ UfX are obtained by a least squares fit. We take the 
N^jy as our estimates of N^ v at the origin of the RNC. Similarly the could be used as 
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an approximation to at the origin of the RNC (though we will not need to do so). The 
appropriate least squares sum is chosen to be 

s(K> ^v) = E E ((U -K- \x\v>w + 

A 4 ij ' 

Note that the symmetry condition AV„ = N\ UfI must be explicitly imposed in this sum. 

Do we have enough data to compute the nine numbers and N^ v l The answer is yes - 
each leg gives us three samples for jVj^ and each computational cell is certain to have more 
than three legs. 

The computation of the K^ v \ a can be performed in a similar manner. The first thing to do 
is to compute the at the origin of each RNC for each computational cell. Now consider 
one specific computational cell. The vertices on the surface of the cell will themselves be the 
origins of neighbouring cells. Each of these cells carries their own estimates of the K^ v in 
their own coordinates. As there is a well defined coordinate transformation between these 
neighbouring frames it is possible to obtain estimates of the K^ v on each of the vertices in 
our chosen cell and in the coordinates of that chosen cell (as indicated in Fig (2)). We can 
then form a linear approximation to each of the six K^ v in the chosen cell. Thus we put 

K nv{x ) = K + K ^ v \ a x 

with the constants K ^ v and K ^ v \ a obtained by another application of the least squares 
method. Finally we take the K^ a as our estimate of K^ a at the origin of our chosen cell. 
Note that in this case we would solve six separate least squares problems, one for each of 
the six K^ v (x a ). In each case we need to compute just four numbers, one K^ v and three 
K^ a . Once again we see that the least squares method is appropriate. 

3. Smooth lattice dynamics 

In the standard ADM 3+1 formulation of vacuum General Relativity there are four constraint 
equations 

= R + K 2 -K^K^ (3.1) 
= % - (3.2) 



and six evolution equations 



00 



dt 



-2NK^ + (3^ + A,,, 
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OK 



a 



where K = K%, = g a/3 R a ^ and R = g^R^. 

Let us concentrate, for the moment, on the evolution equations. When cast as evolution 
equations along the timelike normals these equations take on a particularly simple form, 
namely, 

= ~2NK, V (3.3) 

dK, 



= -^V + N + KK ^ - 2K ^K) (3.4) 

The easiest road to these equations is to simply put the shift vector to zero in the original 
equations. However these equations can also be easily derived by direct calculation. 

The set of equations from (1.2) through to (3.4) can be viewed as a large set of coupled 
ordinary differential equations controlling the future evolution of the lattice. These equations 
will need to be solved by some numerical integration scheme to generate the lattice at 
successive time steps. For pure simplicity we shall demonstrate one time step using the naive 
Euler method, namely, that each quantity is updated according to x <— x + xSt. Note that 
we do not advocate the use of this naive integration scheme in any realistic implementation 
(suitable schemes would include predictor-corrector, leapfrog or Runge-Kutta methods). 

We will start with the lattice data on So, namely, all of the Ly , 0i, fa, dL\j / dt, d6i/dt and 
dfa/dt. We will also assume that well defined rules are given for choosing the N, and 7 M 
on each vertex of the lattice for all time. We will finish with this data fully evolved to the 
next hyper surf ace E^. 

Here is our proposed algorithm. 

1. For each computational cell 

1.1 Use equations (1.2,1.3) to compute the xf,dxf/dt and the R^ a p- 

1.2 Use equations (1.4) and its time derivatives to compute the ccy and datij/dt. 

1.3 Use equations (2.3,2.6) to compute the dL'^/dt and da'^/dt 

1.4 Use equations (2.9,2.10) to compute the 

1.5 Use equations (2.11,2.12) to compute N^ v . 

1.6 Update Lij using dLij/dt. 

1.7 Update xf using dxf/dt. 

1.8 Update 9i,fa using dOi/dt, dfa/dt. 
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1.9 Update K^ v using equations (3.4). 

2. Interpolate back to the primary lattice. 

3. For each computational cell 

3.1 Update dL'^/dt using equations (2.9). 

3.2 Update dLij/dt using equations (2.3). 

3.3 Update da'^/dt using equations (2.10). 

3.4 Update daij/dt using equations (2.6). 

3.5 Update d9i/dt,d(f>i/dt using the time derivative of equations (1.4). 

The are a number of points that need to be made. As each leg will appear in one or 
more computational cells there is an ambiguity in how L y - and dLij/dt should be updated. 
For example, in the cubic lattice of Fig (5), each leg is contained in exactly two different 
computational cells. In this case the L y - and dLij/dt could be updated as the average 
from both cells. Other options are possible and should be explored by direct numerical 
experimentation. Notice that there is no such ambiguity in updating the angles a y - since 
these are defined at the origin of each computational cell. 

The calculations required in step 2 are non-trivial and arise because the shadow lattice will 
drift from the primary lattice (see Fig (3)). As with the computation of the K^ a we would 
need to use coordinate transformations between neighbouring cells and a least squares linear 
approximation to interpolate the K^ v from the shadow lattice back to the primary lattice. If 
one chooses the drift vector to be zero at the origin of each cell then this step is not needed 
and steps 1 and 3 can be merged as one larger step. This is a considerable computational 
advantage. But does it restrict the class of spacetimes that can be built by this method? 
Setting the drift vector to zero everywhere on a lattice is much the same as setting the shift 
vector to zero everywhere on a finite difference grid. We know that in the later case the 
lapse function may need to be carefully chosen so as to avoid the numerical grid collapsing 
on itself (eg. as occurs for iV = 1,N^ = in a Schwarzschild spacetime). The same should 
apply to a lattice - the lapse will need to be carefully chosen to ensure that the lattice does 
not collapse on itself. There are also good arguments for using a shift vector to encourage 
the grid points to maintain a reasonable structure (eg. that the density of grid points does 
not drift too far from some preferred arrangement). It seems reasonable to expect that many 
interesting spacetimes can be explored using a zero drift vector and an appropriately chosen 
lapse function. Despite this, there are known cases (in finite difference calculations) where 
a non-zero shift vector is highly desirable. For example, in the apparent horizon boundary 
condition used by Anninos et al [3] a non-zero shift vector is explicitly used to ensure that 
grid points do not fall into a black hole and to retain a good resolution of the metric in the 
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vicinity of the apparent horizon. Such a condition led to a dramatic improvement in the 
long term stability of the numerical evolution. We can expect that a similar condition will 
be needed in the smooth lattice approach (though this has yet to be tested). 

Another important issue is to what extent does the extensive use of least squares approxi- 
mations (which are required in steps 1.1, 1.4, 1.5, 2, 3.5) effect the accuracy and stability 
of the evolution? The least squares approximations are bound to introduce some degree of 
smoothing in the estimates of N\^ V1 K^ v etc. Does this wash out unimportant numerical 
noise or important short scale variations in the lattice data? Most likely both important and 
unimportant information is lost. The question is does this effect the long term stability of 
the evolution? Again this can be investigated by way of a number of numerical experiments 
on a series of successively refined lattices. 

The reason that we advocated the use of least squares was that it overcame the problems 
of estimating various quantities from an overdetermined system of equations. One could 
imagine a highly idealised case of an overdetermined system which admits a well defined 
consistent solution. For example four points in a plane constitute an overdetermined system 
for a straight line. But properly chosen those four points will determine one straight line. 
In all other cases a straight line approximation could be constructed by a least squares 
method. Suppose now that we wished to impose some dynamics on this toy problem. Let 
us suppose that we started with four points that lie on one straight line (ie. a consistent yet 
overdetermined system). If we impose individual evolution equations for each of the four 
points then those four points may no longer consistently determine a unique straight line. If 
however we impose the evolution equations on the parameters of the straight line then we 
can consistently evolve all four points. This same observation applies to our evolution of the 
smooth lattice. It would be wrong to search directly for evolution equations for each of the 
leg lengths. In our scheme we impose the dynamics on the K^ v which we extract from an 
overdetermined system of equations on the lattice data. Our hope is that if we start with 
a lattice for which we have a good least squares fit for etc. then the evolution will 
continue to yield a good fit. Again this is speculative and must be tested by direct numerical 
computation. 

The discussion so far has dealt entirely with the evolution of the lattice. Some words are thus 
in order regarding the construction of initial data on the lattice. We have already seen how 
all of the quantities needed for the constraint equations can be extracted from the lattice. 
Thus we are able to evaluate the right hand side of the constraint equations. Should this be 
non-zero then clearly we do not have valid initial data and some corrections must be made. 
Following the philosophy espoused in the previous paragraph the corrections to the lattice 
data should be driven by corrections in the coordinate data. How this might be achieved 
is far from clear. Can an approach similar to York's conformal method be used? Part of 
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that method is to propose a 3- metric of the form g^ v = <p g ^ where g ^ is a given 3- metric 
(commonly chosen to be flat) and is the conformal factor. With two metrics one has a 
choice as to which should be expressed in Riemann normal form. If we choose g then we 
can imagine also having a conformal lattice with leg lengths Ly. The Ly of the physical 
lattice could then be estimated as = <pi<f)jLij (as suggested by Piran and Williams [4] ). 
It is unclear at present how York's transverse traceless tensor can be represented on the 
conformal lattice. On the other hand if one chooses to express the physical metric g^ v in 
Riemann normal form then what form does a flat g ^ v take in these coordinates? We conclude 
that it may not be a simple task to apply York's method to a smooth lattice. 

Another important question concerns the order of the discretisation errors. Previous calcu- 
lations (see [1] ) showed that the smooth lattice method yielded 0(e 2 ) accurate estimates 
for both the metric and the Riemann tensors. However, this is unlikely to be the case for 
generic lattices. The formal truncation error in the Taylor series expansion of the metric, 
equation (1.1), is 0(e 3 ) and thus we can expect the discretisation error in the curvatures to 
be of 0{e 1 ). This will couple to the smooth lattice equations leading to an expected error in 
the metric of 0{e 1 ). This is less than optimal - we would prefer an 0(e 2 ) error for generic 
lattices. The challenge then is to find ways in which second order accurate estimates for the 
curvatures on the lattice can be obtained. This also applies to other term such as N^ v and 
Ka/3\v We hope to report on this in a later paper. 

4. An example : The Kasner cosmology 

The Kasner metric 

ds 2 = -{dt'f + f 2 x {t'){dx') 2 + f 2 (t')(dyf + f 2 z (t')(dz') 2 

describes a homogeneous cosmology with a T 3 topology. This metric, for the particular 
choice f x (t') = t Px , f y (t') = t Py , f z (t') = t Pz , is a solution of the vacuum Einstein equations 
when 

1 = Px + Py + Pz = pl + V\ + Pi 

A trivial transformation of the spatial coordinates 

X = fx(t')(x' ~ x' ) 

y = fy(t'W - y'o) 

z = f z (t')(z' - z' ) 

where (x,y,z)' is any nominated point, leads to a 3-metric which, at the chosen point, is 
in Riemann normal form with g^ = diag(l, 1, 1). This fact and the clear simplicity of the 
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metric suggests that the 3+1 smooth lattice approach should yield accurate approximations 
to the Kasner metric. 



In fact this problem is so simple that we will be able to show that there is no discretisation 
error. Of course this fortuitous result can not be expected to occur in other less symmetric 
spacetimes. 

The 3-dimensional lattice is chosen to be a cubic lattice with the addition of two extra 
diagonal legs to each 2-dimensional face. The typical RNC cell will then consist of the 
central vertex, its six nearest neighbours, the six legs that join them and the twelve diagonal 
legs not attached to the central vertex (see Figs (4-5)). The data for each RNC cell will be 
the 6 radial leg lengths L ;, 12 diagonal leg lengths Ly, the 6 angles a y - and their respective 
time derivatives, dL i/dt,dLij/dt and daij/dt. 

Throughout the evolution we will impose the following gauge freedoms. 

o The lapse function has the value one everywhere, Ni — 1. 

o The shift vector is zero at the origin of each RNC, = Ntf. 

o The drift vector is zero everywhere, = 7^. 
While on the initial t = to hypersurface we will assume 

o The diagonal legs are constrained by L 2 - = L 2 oi + L 2 -. 

o The angles (9i,<f>i) are constrained such that Oy = n/2 and daij/dt = 0. 

The choice ay = -ir/2 corresponds to placing the vertices on the coordinate axes while 
da^/dt = ties them to the axes (though they are free to move subsequently off the axes). 
The constraint on the diagonals ensures that the 3-lattice is initially flat, ie. R^ a f3 — 0. 

We will defer stating the homogeneity conditions until we have explicitly calculated all of 
the K^ v . 

Starting with equations (1.2-2.1) we quickly obtain 




(xJ f ) = (L 1+ ,0,0) 
(xJ f ) = (0,L 2+ ,0) 
( a £ f ) = (0,0,L 3+ ) 



(^_) = (-L 1 -,0,0) 
(xJ_) = (0,-L 2 -,0) 
(^_) = (0,0,-L 3 _) 



(A^ + ) = (L 1+ ,0,0) 
(A^ + ) = (0,L 2+ ,0) 
(A^ + ) = (0,0,L 3+ ) 



(7Vf_) = (-V>°>°) 
(iVp = (0,-L 2 _,0) 

= (0,0,-L 3 -) 
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where a dot denotes d/dt. 

Since 7^ = for t > to while dotij/dt = on t = to we find from (2.2-2.3) that 

for t > to 
for t > to 



while 



dx't 


_dx? 


dt 


dt 


dL' tj 


dLij 


dt 


dt 


Mi 


da^ 


dt 


dt 



for t — to 



Using the above we find that equations (2.9-2.10) may be reduced to 

Lij^f = -K^Ax'j (2.90 
= -2K^x^x v j (2.10') 

on t — to- From (2.10') we see immediately that K^ v = for /i 7^ v. The remaining equations 
are of the form 

dL oi ± 2 

°* ~dF ~ ~ KnL oi ± 
dLj±j± 2 2 

L i ± 3 ± ~dT~ = ~ KiiL oi± ~ K 3i L oj± 

for i j drawn from the set {1,2,3}. We can see that the second equation is a trivial 
consequence of the first equation and the assumed constraint on the diagonals. Thus this 
second equation is redundant and may be excluded when solving for the K^ v on t = to- Thus 
we obtain the following equations for the K^ v 

dL Ql - 2 

L oi-^r-' KnL oi- 

dL - 2 

L M-— = - K22L 02- 

dL nn - 



L 0l+ 


dL Q1 + 


= 


dt 


L 02+ 


dL 02+ 


= -K 2 2L 2 m+ 


dt 


L 03+ 


dL Q3+ 




dt 



This is clearly an overdetermined system for Kn,K22 and -K33. One could choose to take 
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an average solution, such as, 



K22 = -r 



1 dL ru + 1 dL ru - 



K33 



L 01 + 


dt 


+ 


L 01~ 


dt 


1 


dL Q2+ 


+ 


1 


dL Q2 - 


L 02 + 


dt 


L 02~ 


dt 


1 


dL Q3+ 


+ 


1 


dL QS - 



, 03 + dt L Q3 - dt 



(which we would also get from a least square or singular value decomposition). However, we 
have yet to impose precisely the condition that the lattice be homogeneous. 

Consider two neighbouring vertices p and p and their respective RNC cells. Let x^ and x^ be 
the RNC coordinates for the vertices p and p respectively. We can align the coordinates such 
that the two vertices lie on, say, the x-axis. We can now impose homogeneity by demanding 
that g^ v and K^ v of the RNC for p be obtained by Lie dragging the and K^ v of the 
RNC for p along the integral curves of d/dx. For our lattice the 3-metric is flat everywhere 
so this condition reduces to just a statement about the extrinsic curvatures, 

This now forces us to make an alternative choice in solving the above system for K^, namely, 

— 1 g?L q1 + —1 dL Q1 - 
Lq^+ dt ^oi~ dt 
^ — 1 dL^2+ — 1 dL^ 2 - ^ ^ 

K33 = ^- = ^- (4.3) 

Lq^+ dt -Zvqj— dt 

The initial data must be chosen to obey these constraints. 

A further consequence of this definition of homogeneity is that 

= K^ j0 , = K^ v \ a 

the second equality following from the RNC condition that = at the origin. Thus we 
see that the momentum constraint is identically satisfied. Since K^ v is diagonal we also find 
that the Hamiltonian constraint (3.1) is just 

= ^11^22 + ^11^33 + ^22^33 (4.4) 
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The main evolution equations (3.4) are now 



dKi 



l 



dt 

dK22 

dt 
dK 33 
dt 

while 



(-K u + K 22 + K 33 )K U (4.5) 
(+K U - K 22 + K 33 )K 22 (4.6) 
(+K U + K 22 - K 33 )K 33 (4.7) 



^ = (4.8) 
dt 

for n^v. 

We can see that these equations, when coupled with (4.1-4.3), provide second order evolution 
equations for the L ;. We also need second order evolution equations for a y - and L y -. 

Using K^y = 0, dK^jdt = on t = to for \i ^ v we find, by differentiation of (2.10'), that 

d 2 aj 



dt 2 

on t = to, while differentiation of (2.9') leads to 



= (4.9) 



also on t = to- These are the required evolution equations for and Ly. These equations 
are nothing more than the second time derivatives of the original constraints on the angles 
and diagonals, namely, ay = tt/2 and L 2 - = L 2 oi + L 2 -. 



We are now in a position to fully analyse the future evolution of the lattice. For our lattice 
the initial data are L j, L y -, cty and their first time derivatives. The complete set of evolution 
equations are (4.1-4.3,4.5-4.10). There is only one initial value constraint, equation (4.4). 

From equations (4.9) and (4.10) we see that our initial constraints on ay and the diagonals 
Lij are preserved by the evolution equations. Equation (4.8) shows that K^ u remain diagonal. 
The dL i/dt are updated from equations (4.1-4.3) and thus the lattice remains homogeneous. 

In summary, the smooth lattice equations guarantee that the future 3-dimensional lattices 
remain flat and homogeneous. 
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By analogy with the Kasner form of the metric we might seek a solution of the form 



L 01+ 


= L 01~ 


= At Px 


L 02 + 


= L 02~ 


= Bt p y 


L 03 + 


= L 03~ 


= Ct Pz 



where A, B, C are constants. It can seen that this is a solution of the constraint (4.4) and 
the evolution equations (4.5-4.7) provided 

1 = Px + Py + Pz = Pi + p\ + p\ 

in exact agreement with the Kasner metric. 

We conclude that for this simple model there is no discretisation error. It should be noted 
that the same observation would occur if one had employed a finite difference approach. This 
is a rather trivial example as it does not touch upon some of the more delicate aspects of 
the method. The extrinsic curvatures where computed without resort to least squares and 
the assumption of homogeneity allowed us to set N^ v and K^ v \ a to zero. There was also no 
need to interpolate the future values of K^ v back to the updated lattice. In a more general 
setting each of these aspects can be expected to require some attention. 

5. Discussion 

There have been other attempts to adapt the ADM 3+1 equations to a lattice. In particular 
there are the works of Piran and Williams [4] and Friedmann and Jack [5] . In both cases the 
equations were presented in the context of the Regge Calculus. This is another lattice method 
in which the metric is assumed to be flat inside each pair of adjacent tetrahedra. The principle 
difference between their approach and ours is in the way in which the ADM equations were 
imposed on the lattice. In Piran+Williams and Friedmann+Jack the equations of motion for 
the lattice were obtained from the ADM 3+1 action principle. They were unable to evaluate 
the action directly because in the Regge Calculus one does not have direct access to quantities 
such as R^p^N^ etc. Instead one has pure scalar quantities such as the defect angle, the 
areas and volumes of simplices. Thus in developing their equations both pairs of authors 
needed to make many (reasonable) assumptions about how various terms in the standard 
ADM action could be translated into a Regge form. This is a somewhat adhoc method and 
thus casts some doubt on the validity of the final equations. Indeed the equations given by 
Friedmann and Jack differ from those given by Piran and Williams (though they do agree 
for appropriate choices of lapse and shift). 
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In contrast our equations have been derived in a systematic manner. The use of Riemann 
normal coordinates has allowed us to extract in a very natural way all of the relevant coor- 
dinate data from the lattice data. This in turn has lead to a very natural adaptation of the 
ADM 3+1 equations to the lattice. Despite this there are still many variations that need to 
be explored. Some of these are of a minor numerical nature (eg. should all the dLij/dt be 
used in estimating K^ v or just a limited set?). Others are potentially much more important, 
in particular does the use of overdetermined systems lead to errors and instabilities in the 
evolution of the lattice? There is also the major question of how one can solve the initial 
value constraint equations. Despite these concerns, we believe that the basics of our algo- 
rithm are well founded and should survive in any subsequent lattice method for numerical 
relativity. 

However, the proof of the pudding is always in the eating. The method has been shown to 
have no discretisation error for the rather trivial example of the Kasner cosmology. This is 
encouraging but it must be noted that this example avoided many of the potential problems 
associated with this method. So it is imperative that the proposed method be put to a 
serious test. We have already successfully constructed the initial data for a time symmetric 
initial data slice in the Schwarzschild spacetime [1] . We are currently applying our proposed 
method to evolve that data. The results will be reported (soon) in a later paper. 
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Figure 1. This figures displays the drift vector 7^ and the relationships between 
the coordinates on the primary and shadow lattices. Clearly 5xf = (— 7^ + Nf*)5t 
leading directly to equation (2.1). We also have, by construction, that xf = x'l 1 on 
Sq. Thus 5xf = dx'- 1 + rfSt which in turn leads to equation (2.2). 



Figure 2. In the first pass over the lattice, K^ v is computed at each vertex (•). In 
a second pass, K^ a at the central vertex is estimated by fitting the linear approx- 
imation K^ v (x) = K^ v + K fIU \ a x a to the data on each vertex. This will require a 
coordinate transformation from neighbouring cells to get data on the boundary of 
the cell. The N^ v at the central vertex can be estimated by first forming estimates 
for iV| M at the centre of each leg (■). That data is then approximated by a linear 
function N^(x) = + N\^x v . We can then use N^ v as an estimate of at 
the central vertex. 



Figure 3. After one time step the shadow lattice has drifted relative to the primary 
lattice. However the updated are recorded on the shadow lattice. Thus an 
interpolation from the shadow lattice back to the primary lattice is required. This 
can be done by forming a linear approximation K^ v (x a ) = K^ v + K^ a x a . and 
evaluating it at x a = 0. Note that no such interpolation is required if the drift 
vectors are set to zero. 



i z 




Figure 4. A typical computational cell in the cubic lattice. This cell consists of eight 
tetrahedra attached to the central vertex. The full cubic lattice can be obtained by 
replicating this structure throughout the space. Note that this leads to overlapping 
tetrahedra. This typical cell can be used in many spaces other than just the Kasner 
cosmology considered in the text. The T 3 topology of the Kasner cosmology can 
be obtained by identifying opposite ends of the cubic lattice. Note that, for clarity, 
some of the legs have not be shown in this figure (eg. (l + 3~), (1~2 + ) etc.) 



Figure 5. An rry-cross-section of (part of) the cubic lattice. The heavy solid lines 
denote the legs of one computational cell. Notice that the diagonal legs attached 
to the central vertex are not included in the computational cell. Note also how 
neighbouring cells overlap each other. 



